Effective single-band models for strongly interacting fermions in an optical lattice 
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To test effective Hamiltonians for strongly interacting fermions in an optical lattice, we numerically 
find the energy spectrum for two fermions interacting across a Feshbach resonance in a double well 
potential. From the spectrum, we determine the range of detunings for which the system can be 
described by an effective lattice model, and how the model parameters are related to the experimental 
parameters. We find that for a range of strong interactions the system is well described by an effective 
t — J model, and the effective superexchange term, J, can be smoothly tuned through zero on either 
side of unitarity. Right at and around unitarity, an effective one-band general Hubbard model is 
appropriate, with a finite and small on-site energy, due to a lattice-induced anharmonic coupling 
between atoms at the scattering threshold and a weakly bound Feshbach molecule in an excited 
center of mass state. 



I. INTRODUCTION 

Systems of ultracold fermionic atoms interacting in an 
optical lattice potential via a Feshbach resonance provide 
remarkable opportunities to realize a zoo of lattice Hamil- 
tonians in a clean and controllable fashion [3-[3]- This 
allows a new tool to study many lattice models famil- 
iar from condensed matter physics such as the Hubbard, 
t — J, and XXZ models, to name a few. In addition to 
the prospect of studying these paradigmatic models ex- 
perimentally in the absence of unwanted complications, 
ultracold gases provide the exciting possibility of gain- 
ing new insight into the proper description of strongly 
correlated systems. 

While it is well-known that a weakly interacting gas in 
an optical lattice can be described by the one-band Hub- 
bard model [j], IH , the situation for strongly interacting 
gas near a Feshbach resonance is much more complicated. 
In the strongly interacting regime, the conventional as- 
sumptions for derivation of the one-band Hubbard model 
for two-component fermions obviously no longer applies, 
since the on-site interaction energy becomes greater than 
the bandgap of the lattice and the off-site interaction 
gets comparable with the atomic tunneling rate @, @- 
However, Refs. @, provide general arguments to show 
that in this case we can still derive an effective single- 
band lattice Hamiltonian, in the form of either a general 
Hubbard model with possibly particle correlated hopping 
rates, or a t — J or XXZ model under different limiting 
situations. This simplification relies on several observa- 
tions: first, with the assumption that the average num- 
ber of atoms per lattice site n < 2, it is unlikely for 
more than two atoms to occupy the same site because 
they are energetically unfavorable due to the Pauli ex- 
clusion and the strong on-site interaction Q. Second, 
with two atoms on the same site, although they will 
populate many lattice bands due to strong interaction, 
the two-atom eigen-levels (called the dressed molecules) 
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have large energy splitting between them 0, H, H4U , and 
for low-temperature physics, only one of these dressed 
molecule levels will be relevant (the other levels can be 
adiabatically eliminated in the derivation). So the effec- 
tive Hilbert space is severely restricted: each site may be 
empty, or populated with either one atom, or one dressed 
molecule with a fixed internal state. Based on restriction 
of the effective Hilbert space and the SU(2) symmetry for 
the underlying physical process, it is derived in Ref. @ 
that the effective lattice Hamiltonian takes the form of a 
general Hubbard model with possibly particle correlated 
hopping rates. The multiple-band population and the off- 
site interaction are taken into account through renormal- 
ization of the effective parameters in the final single-band 
model. In different interaction regions for the atoms (con- 
trolled by the external magnetic field through the Fesh- 
bach resonance), the dimension of the effective Hilbert 
space on each lattice site can be further reduced, leading 
to an effective t — J model for the atoms or an XXZ 
model for the dressed molecules [2,6]. Note though, that 
in these regions the t — J or the XXZ models in general 
do not arise from the Hubbard (or the general Hubbard) 
model via perturbation theory, and they should be taken 
as the basic models for the underlying physical process. 

This paper serves two purposes: First, we provide nu- 
merical tests for the effective lattice Hamiltonians ap- 
plied to the two-body problem. The effective Hamilto- 
nian should describe two-body physics as well as many- 
body physics. For the two-body problem, we can solve it 
exactly in a double well lattice from the basic strongly in- 
teracting field Hamiltonian and calculate the low-energy 
eigen-levels of the system. We then compare the results 
from the exact numerical solutions with prediction from 
the effective single-band models. This comparison deter- 
mines the validity region for each effective single-band 
model. We note that the two-body strongly interact- 
ing problem has also been solved for a three-dimensional 
lattice in the BCS region [12j, for a one-dimensional lat- 
tice (l3| . and for two bosons in a triple- well potential 
14], based on different kinds of numerical methods. Sec- 
ond, with comparison of the two-body physics, we can 
microscopically derive the parameters for the effective 
single-band models from the relevant experimental pa- 
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rameters. The determination of the effective model pa- 
rameters is very important for application of these mod- 
els to do quantitative calculations to compare theory and 
experiments. We find that for a wide range of strong in- 
teraction, the system is well described by an effective 
t — J model, and the effective superexchange term, J, 
can be smoothly tuned through zero on either side of 
unitarity. Right at and around unitarity, an effective 
one-band general Hubbard model is appropriate, with a 
finite and small effective on-site energy, due to a lattice- 
induced anharmonic coupling between atoms at the scat- 
tering threshold and a weakly bound Feshbach molecule 
in an excited center of mass state (the so-called anhar- 
monic induced resonances, see Ref. fl5j). 

This paper is arranged as follows: In Sec. II we spec- 
ify our numerical approach for exact solution of the two- 
body problem in a double-well potential, which is based 
on the stochastic variational method [16}, a widely-used 
method in nuclear physics for solution of the few-body 
problems. The main results of this paper are shown in 
Sec. III. First, we specify the validity regions for each 
effective single-band lattice model across the whole in- 
teraction region, from the BEC limit, to resonance, and 
to the BCS limit. Second, we derive the renormalized 
parameters for these effective lattice models. The values 
for some of these parameters are pretty counter-intuitive, 
and in such cases we discuss their physical origin and con- 
sequence. 



II. METHODS 

We consider two distinguishable fermions of mass m in 
an external potential V (x) , interacting via a short range 
potential U (r) characterized by its s-wave scattering 
length, a s . In experiments, the external potential is typ- 
ically the sum of a periodic potential Vq JliLi cos2 k^Xi 
and a harmonic confining potential. To make our calcu- 
lations easier we will model this by considering a double- 
well potential along the z-axis (formed by Taylor ex- 
panding cos 2 ki,z) and a harmonic potential in the other 
two transverse directions, with the frequency, uj, cho- 
sen such that the potential is locally isotropic at the 
bottom of each well, as shown in Fig. [TJ In fact, this 
double- well problem is interesting in its own right, as it 
is relevant to experiments with gases confined in optical 
superlattices [17||. Since the barrier is the same as for 
the full lattice, the hopping rates should be nearly un- 
affected. Also, the on-site interaction energy should not 
be affected qualitatively. The only relevant information 
that the double-well approximation renders inaccessible 
is rates for any next-nearest neighbor process. The effec- 
tive single-band lattice Hamiltonians we are comparing 
belong to the tight-binding models with interaction only 
among the nearest neighbors, so it is enough to consider 
a double-well potential to test these Hamiltonians and 
their parameters. 

Due to the harmonicity of the transverse trap, the cen- 




FIG. 1: (a) The double- well potential along z modeling a 
periodic potential; (b) contour plot of the locally isotropic 3D 
double-well potential. 

ter of mass (cm.) motion in the transverse direction sep- 
arates out and is thus neglected in the rest of the dis- 
cussion. However, along the axis of the double-well, the 
cm. motion is not separable from the relative motion. 
The two atom system then has three relevant coordi- 
nates: the relative coordinates z and p, along the axial 
and transverse directions, respectively, and Z ', the axial 
cm. coordinate. In terms of these coordinates, the ex- 
ternal potential approximating an optical lattice of depth 
V is V (p, z, Z) = V (p) + V (z, Z), where 

V(p) = V k 2 L p 2 /2 

viz z) Vn V HD" 1 ^ 1 - 2 ") !: 2n ( z , ^ 2n 
v(z, z) -Vo}^ T{1 _ 4n) r{1 + 4n) k L [Z±-j 

(i) 

and r (x) is the Euler gamma function. 

The exact form of the interaction, U {^J p 2 + z 2 ^j , is 

irrelevant in the low-energy limit as long as its range 
is much smaller than the average interatomic distance 
and the trap length scale, and most analytical treat- 
ments use a zero-range pseudopotential. Numerically, it 
is easiest to use a finite-range attractive Gaussian interac- 
tion U (r) — — [/oexp (— r 2 /r 2 ), where we typically take 

r = O.Qby/h/mu}. Finite-range effects should be neg- 
ligible for such small values of ro, and we have verified 
this by repeating the calculations with ro = O.ly/h/mu}. 
The free space scattering length is varied by adjusting 
the strength of the interaction, Uo. 

Adopting units such that ki, = 1 and En — 
ti 2 k 2 j /2m = 1, the Hamiltonian may be written as 

2 d d n d 2 Id 
p dp dp dz 2 2 dZ 2 

+ 2 1 4 + V(p,z,Z)-U e-( z2+ P 2 )^ (2) 
P 

where mg is the relative angular momentum, which is a 
good quantum number due to axial symmetry. In the 
following we will only consider mi = 0, since, in the limit 
as r goes to zero, the interaction does not affect states 
with m( ^ 0. 
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We find the low-lying states of the system using a 
stochastic variational method 11611 recently introduced to 
the ultracold gas community [13]. In this approach, the 
variational wavefunction takes the form 

N 

* (p, z,Z) = J2 <*i exp (~p 2 /4 - z 2 /bl - Z 2 / C 2) , (3) 

i 

where a is a linear variational parameter, {a, b, c} are 
nonlinear variational parameters which define the basis 
elements, and N is the size of the basis set. The non- 
linear parameters are selected from stochastically gen- 
erated pools of candidates to minimize the variational 
energy {9 a bc\H \^ a bc) / abc\^ abc) ■ The basic algorithm 
is as follows: starting with a set of N — 1 basis states, 

1) a pool of (in our calculations) 25 new basis states is 

randomly generated, each defined by a given value 
of {a^b^a}; 

2) for each of the 25 possible TV-dimensional basis sets 

formed by adding one basis from the candidate 
pool, the energy is minimized with respect to a; 

3) the new basis set that yields the lowest energy is kept 

and the previous steps are repeated until the basis 
size, N, increases to the desired number. 

Once every few iterations, the existing basis set is op- 
timized by the following refining process: starting with a 
set of N basis states and n = 1, 

A) a pool of 25 replacement basis states is randomly gen- 

erated, each defined by a given value of {a n , b n , c„}; 

B) for each of the 25 possible iV-dimensional basis sets 

formed by replacing the n th old basis state with 
a new one from the candidate pool, the energy is 
minimized with respect to a; 

C) if the lowest of these 25 energies is lower than the 

current variational energy, the n th old basis state 
is replaced by the new optimal one and the previous 
steps are repeated for n — 1...N. 

The great advantage of this stochastic variational 
method is that it is able to avoid getting stuck in the local 
minima of the energy landscape that plague determinis- 
tic variational methods. Unlike Monte Carlo algorithms, 
though, there is no well-characterized statistical error. 
The resulting variational energies set upper bounds on 
the true eigenenergies, but we cannot attach a rigorous 
error estimate; we only check that the results seem well- 
converged. For a detailed discussion of the method, see 
Ref. [II]. 

We typically achieved fairly good convergence for N ~ 
300. Although in principle the nonlinear basis optimiza- 
tion must be performed for each value of a s , actually 
the basis set does not change too much as one sweeps 
across resonance except to include narrower and narrower 



Gaussians for positive a s where deeply bound molecules 
form. Apart from deeply bound states, the change in 
the wavefunction is mainly due to changing the expan- 
sion coefficients, a. To save computational time then, 
we performed the nonlinear basis optimization for four 
different values of a s across resonance, joined the four 
optimized basis sets, and simply minimized the energy 
with respect to a using the resultant basis set of about 
1200 elements for all values of a s . As a result, very deeply 
bound energy levels may not be fully converged, but we 
are only interested in the energy range around the lowest 
non-interacting levels. In this range, our results appear 
to be converged. 



III. RESULTS 

A. Regions of model validity 

In Figs. [5] and [31 we show the spectrum of two atoms 
interacting in a double-well potential across a Feshbach 
resonance for Vq = 8E r and Vq — 10E r . For clarity, we 
have diabatized the spectrum for — l/fc^a s < —4, omit- 
ting the plunging levels. We have also omitted the ex- 
actly flat, non-interacting levels. In the absence of trap 
anharmonicity, there are three kinds of curves present: 
plunging levels corresponding to tightly bound molecules 
in motionally excited states, flat levels corresponding to 
atoms in separate wells, and sigmoidal levels correspond- 
ing to interacting extended atom pairs. Levels corre- 
sponding to states of similar parity never cross, instead 
undergoing a rich set of narrow avoided crossings. These 
are due to Feshbach-type resonances induced by the an- 
harmonic coupling of the center of mass and relative mo- 
tion (lsj . Away from unitarity, these avoided crossings 
become even narrower and are irrelevant. Some of our 
calculations of highly excited even molecule states are 
apparently not as well converged as the corresponding 
calculations for the odd states, since the even states (solid 
lines) should be the lowest of each plunging doublet as 
the molecules become tightly bound and act as a single 
particle. However, this is not important for our purposes. 

To approximate the two-atom, double-well physics 
with a general Hubbard model requires a separation of 
energy scales such that there is a manifold of four energy 
levels well-separated from all the others, corresponding 
to a doublet of doubly occupied states and a doublet 
of singly occupied states. This condition is satisfied far 
from unitarity on either side. (On the scale of the plots, 
it may be hard to distinguish the two levels in the doubly 
occupied doublet, as their splitting is on the order of the 
dimer tunneling energy. The levels in the singly occupied 
doublet are likewise very close, with splitting on the order 
of the superexchange energy, but we have only shown the 
energy of the singlet state since the triplet state is odd 
in the relative coordinate and thus independent of the 
interaction.) We have marked the low-energy regions of 
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FIG. 2: (Color online.) (a) Spectrum of two interacting atoms 
in a three-dimensional double-well potential vs. inverse free 
space scattering length with Vo = 8Er. Solid (dashed) lines 
correspond to states of even (odd) symmetry in Z. Only 
the first few plunging levels are shown, (b) Close-up of the 
strongly interacting region. The noninteracting states odd in 
z are not shown. 



model validity using a rule of thumb that the energy sep- 
aration between the relevant low-lying manifold and the 
nearest level outside the manifold should be at least five 
times larger than the energy range of the manifold. Ac- 
tually, in enforcing this requirement, we have taken into 
account that the discrete levels turn into bands of width 
~ At when extending the two-site potential to an infinite 
one-dimensional lattice. 

The small disconnected regions near unitarity where a 
general Hubbard model is applicable, shown in Figs. 2(b) 
and |3(b)| are qualitatively different from the weakly in- 
teracting regions. Here the on-site dimers correspond to 
Feshbach molecules in an excited center of mass band. 
The usual ground band Feshbach molecules are far- 
detuned and irrelevant. The coupling to the relevant 
singly occupied states (which are in the lowest center of 



FIG. 3: (Color online.) Same as Fig. [2] but for Vo = WE R . 



mass band) is facilitated by the anharmonicity of the po- 
tential. This is a very interesting phenomenon whereby 
the presence of a narrow anharmonicity induced reso- 
nance [l5j near the wide free space Feshbach resonance 
allows for an effective single-band description where the 
on-site energy becomes very small instead of arbitrarily 
large. As higher excited molecular bands become rele- 
vant, eventually the molecule bandwidth becomes com- 
parable to the bandgap and one can no longer apply a 
general Hubbard model. This is the case in Fig. |2(b)| at 
-l/k L a s - -1.7. 

Likewise, approximation by a t — J model requires that 
two singly occupied states are well-separated from the 
others, and these regions are also marked in Figs. [2] and 
[3l (Again, one of these states is not shown and would be 
indistinguishable anyway on the scale of the plot.) Note 
that the form of the t — J model is valid at detunings 
where the general Hubbard model is not. There the t — J 
model does not come from the usual perturbative treat- 
ment of the Hubbard model for \U\ ^> t, which gives 
J = 2t 2 /U [l9|. Instead it is the fundamental model in 
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that region and there are no a priori constraints on the 
parameters. The regions of validity shown in Figs. [5] and 
[3] for the two types of lattice models are one of the main 
results of this paper. 

We have performed the calculations for lattice depths 
of V = 2 — IOEr. For Vq < 4E R , a lattice model is not 
valid for any value of the scattering length. At Vq — QEr, 
the bandgap has increased enough that the lattice mod- 
els become valid in regions, similar to those plotted in 
Fig. [2j but narrower. These regions quickly widen when 
the lattice depth is increased to Vb = 8Er. As the lat- 
tice is deepened to Vq — IOEr, the t — J model contin- 
ues to quickly expand its validity further into the weakly 
interacting regime, overlapping with the regions where a 
Hubbard model is valid. It expands only slightly closer to 
unitarity. A general Hubbard model also becomes valid 
over only a slightly larger region. 



B. Model parameters 

Above we have determined the structure of the relevant 
model Hamiltonians across unitarity. Now we consider 
the details of the Hamiltonians and their parameters. Ne- 
glecting off-site interactions, the general Hubbard model 
takes the form m 



Hghm = 53 (rats- A - /-0 n i 

/ 

-E 



t+(g - t) {n is + n js ) 



+ (tda + t-2g) n is n ]S 



(4) 



where the prime on the sum means that only near- 
est neighbor terms are included, a denotes a fermion 
component {t, 4-} an d o denotes the other component, 
riia- = a icr aicr, A is the on-site interaction energy, [i is 
the chemical potential, g is the particle-assisted tunnel- 
ing rate, tda is the rate for an atom to hop from a dou- 
bly occupied site to a singly occupied site, and td is the 
dimer tunneling rate. For just two atoms on two sites, 
the Hamiltonian can be written as the matrix 



H, 



(2,2) 
GUM 



/A —g -g -t d \ 

-g o -g 

-9 -g 

\-*<j -9 -g A 



2fil 



(5) 



in the basis{|tl,0),|t,l),U,t>,|0,n>}. 

From the spectrum we can extract the model param- 
eters as a function of the experimental parameters by 
choosing them such that the lattice model correctly re- 
produces the four energy levels. The resulting values de- 
pend on the scattering length and are shown in Fig. 3] for 
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FIG. 4: (Color online.) (a) Particle-assisted hopping rate and 
(b) on-site interaction energy vs inverse scattering length for 
Vo = 8Er. Dimer hopping is negligible and not shown. 



Vq = 8Er. The model is not valid in the hatched regions, 
in accord with Fig. [2] The particle-assisted hopping rate, 
g, generally differs from the single-particle hopping rate, 
t, by about 10%. (We use the value of t obtained from 
the one- atom spectrum.) Away from resonance, the on- 
site interaction energy, A, is 25% larger than the stan- 
dard Hubbard U computed from the overlap of the lowest 
Wannier functions, but this is no doubt due in part to 
the additional confinement compared to the infinite lat- 
tice case. We have not shown the dimer hopping rate 
away from resonance, as it is negligible - less than 1% 
of the single-particle hopping rate in this case. Far from 
unitarity, where the interaction is too weak to populate 
higher bands, the general Hubbard model reduces to the 
one-band Hubbard model. In general, except for within 
the slivers of model validity around the anharmonicity 
induced resonances [HI , Fig. [4] shows that the general 
Hubbard model is qualitatively similar to the standard 
Hubbard model with some quantitative corrections as one 
begins to approach the strongly interacting region. 
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FIG. 5: (Color online.) (a) Particle-assisted hopping rate, 
dimer hopping rate, and (b) on-site interaction energy vs in- 
verse scattering length in the narrow region around the first 
anharmonicity induced resonance where a general Hubbard 
model is valid for Vb = 8Er. 
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FIG. 6: (Color online.) Same as Fig. [5] but in the narrow 
region of model validity around the second induced resonance. 



FIG. 7: (Color online.) Superexchange energy vs. inverse 
scattering length for Vb = 8Er. 



The effective t — J model has the familiar form 
H t -.j = -/iy^rw 



hi 



t ^Z P a a \a a 3<y P a ~ J (si • Sj - UiTlj / A) 



(6) 



Around the first anharmonicity induced resonance, 
though, the general Hubbard model is strikingly different 
than the single-band Hubbard model, as shown in Fig. [5] 
Of course, the physical single-band approximation is not 
expected to hold in this strongly interacting region, but 
it is still interesting to contrast the two models. Most 
notably, the effective on-site interaction does not become 
unbounded near the free space Feshbach resonance. Al- 
though the relevant on-site dressed molecule state resem- 
bles the lowest dressed molecule state in that it is tightly 
bound internally, it is only slightly detuned energetically 
from the atom pair singlet state due to its excited center 
of mass motion. Also, the dimer tunneling is no longer 
negligible, as might be expected for an excited dimer in 
a relatively weak lattice. Only the particle-assisted hop- 
ping is qualitatively similar to that plotted away from 
unitarity in Fig. HI 

Similarly, for the small region around the crossing of 
the second excited molecule state and the lowest atomic 
state, one can obtain the general Hubbard model param- 
eters shown in Fig. [SI The particle-assisted hopping and 
the on-site interaction energy are much the same as in 
Fig. [SJ However, here the dimer tunneling rate is much 
larger, twice as large as the atom tunneling rate, due to 
the excited nature of the relevant molecule state. 



where P a is a projector onto the subspace with at most 
one atom per site, J is the superexchange energy, m — 
J2a n w and s i = J2a,a> <H<r a W CLia' / 2 , with er denoting 
the Pauli matrices. For just two atoms on two sites, the 
Hamiltonian can be written as the matrix 

HSS = - (jf) - 2,1 (7) 

in the basis {| t)4)> I 4-jt)}- We extract from the spec- 
trum the effective superexchange energy relevant for the 
regions where a,t — J model is valid. We have plotted this 
in Fig. [3 along with the result of 2t 2 /U one would derive 
from the conventional one-band Hubbard model. Again, 
this is just for reference, as there is no expectation that 
the physical single-band approximation is valid in the 
strongly interacting regime. It is interesting to note that 
J does not vanish at unitarity, but in fact passes through 
zero on both sides of unitarity (in fact, more than once 
on the positive-a s side). This implies that one may find 
some transition to an exotic phase there, where the in- 
teraction is dominated by higher-order or next-nearest 
neighbor processes. A similar zero-crossing has previ- 
ously been found on the positive g s side of resonance in 
a completely different calculation [20|. 

As the lattice depth is adjusted, the behavior of the pa- 
rameters remains qualitatively unchanged and the maxi- 
mum values of the particle-assisted hopping rate and the 
superexchange energy remain more or less the same. The 
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only exception is that td/t decreases quickly as the lat- 
tice becomes deeper, as expected. Where both models 
are valid, the t — J model can be derived from the gen- 
eral Hubbard model and J = 2g 2 /A. 

IV. CONCLUSIONS 

We have discussed the basic form of candidate effec- 
tive single-band lattice Hamiltonians to describe the low- 
energy physics of ultracold fermionic atoms in an optical 
lattice. General considerations of the relevant Hilbert 
space and system symmetry Q lead to two possibilities 
for energies near the non-interacting ground state: an 
effective single-band generalized Hubbard model and a 
t — J model. 

We have performed numerical calculations of the spec- 
trum of two interacting fcrmions in a double- well poten- 
tial to determine under what conditions one of these lat- 
tice models is a good description of a physical system 



in a periodic optical potential. By requiring the lattice 
models to reproduce the low-energy two-site two-atom 
physics in their respective regions of validity, we have 
determined the on-site and nearest-neighbor parameters 
of the models. We find that at unitarity there exists a 
valid effective single-band Hubbard model with counter- 
intuitive weak on-site interaction. We also find that near 
unitarity there exists a valid t — J model whose superex- 
change energy can be tuned through zero on either the 
attractive or repulsive side of resonance. These models 
should prove useful starting points for future theoreti- 
cal or experimental considerations of strongly correlated 
many-body physics in an optical lattice. In particular, 
the ability to tune J through zero suggests an interesting 
phase diagram in the vicinity of the crossing. 

This work was supported by the AFOSR through the 
MURI quantum simulation program, the DARPA OLE 
program, the IARPA, and the ARO through a MURI 
program. 
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